import numpy as np
import csv
from datetime import datetime
from time import strftime
import matplotlib.pyplot as plt
import math
import random
from scipy import stats
#import ols
from scipy.stats import f
import statsmodels.api as sm


def getCurTime():
    """
    get current time
    Return value of the date string format(%Y-%m-%d %H:%M:%S)
    """
    format='%Y-%m-%d %H:%M:%S'
    sdate = None
    cdate = datetime.now()
    try:
        sdate = cdate.strftime(format)
    except:
        raise ValueError
    return sdate

def build_data_list(inputCSV, datatype="int"):
    sKey = []
    fn = inputCSV
    f = open(inputCSV)
    #ra = csv.DictReader(file(fn), dialect="excel")
    ra = csv.DictReader(f, dialect="excel")

    if datatype=="int":
        for record in ra:
            #print record[ra.fieldnames[0]], type(record[ra.fieldnames[-1]])
            for item in ra.fieldnames:
                temp = int(float(record[item]))
                sKey.append(temp)
    elif datatype=="float":
        for record in ra:
            for item in ra.fieldnames:
                temp = float(record[item])
                sKey.append(temp)
    else:
        print 'incorrect value in second parameter (datatype)'
        return


    sKey = np.array(sKey)
    sKey.shape=(-1,len(ra.fieldnames))
    return sKey

def geneBpChoices(tmin, tmax, num):
    tmin = float(tmin)
    tmax = float(tmax)
    num = int(num)
    step = (tmax - tmin)/num
    tbpchoices = []
    for i in range(num):
        tbpchoices.append(tmin+i*step)
    return tbpchoices

def pickBreakpoint(response, predictor, bpvarible, bpChoices, predictorName):
    results = np.zeros((len(bpChoices)-1, 2))
    print bpChoices
    bpvarible = np.array(bpvarible)
    predictor = np.array(predictor)
    predictor_t = np.transpose(predictor)
    predictorName = np.append(predictorName, 'bpvarible')
    #print predictor_t, predictor_t.shape
    
    for i in range(len(bpChoices)-1):
        x2star = (bpvarible - bpChoices[i]) * np.greater(bpvarible, bpChoices[i])
        tempPredictor = np.append(predictor_t, x2star)
        tempPredictor.shape = (len(predictorName),-1)
        tempPredictor = np.transpose(tempPredictor)
        #print tempPredictor
        tempmodel = ols.ols(response, tempPredictor,'y', predictorName)
        results[i,0] = i
        results[i,1] = tempmodel.R2
    print results

    optBP = int(results[np.argmax(results, axis = 0)[1],0])
    print 'Optimal Index:', optBP
    print 'Optimal changepoint: ', bpChoices[optBP], ' exp value: ', np.exp(bpChoices[optBP]), ' with R2 = ', results[optBP, 1]

    x2star = (bpvarible - bpChoices[optBP]) * np.greater(bpvarible, bpChoices[optBP])
    tempPredictor = np.append(predictor_t, x2star)
    tempPredictor.shape = (len(predictorName),-1)
    optPredictor = np.transpose(tempPredictor)
    optmodel = ols.ols(response, optPredictor,'y', predictorName)
    #optmodel = ols.ols(response, optPredictor,'y',predictorName)
    
    y_hat = np.dot(optmodel.x, optmodel.b)
    #fileLoc = filepath + 'pieceR_model2_y_hat.csv'
    #np.savetxt(fileLoc, y_hat, delimiter=',', fmt = '%s')
    
    print optmodel.summary()
    print 'MSE =', optmodel.sse
    #print 'before bp, b0 =',optmodel.b[0], ', b1 =', optmodel.b[1], ', bd =', optmodel.b[2]
    #print 'after bp, b0 =', optmodel.b[0] - optmodel.b[3] * bpChoices[optBP], ', b1 =', optmodel.b[1], ', bd =', optmodel.b[2] + optmodel.b[3]

    #calpredictedvalue(zip(x1, predictor), bpChoices[optBP], zip(optmodel.b, [optmodel.b[0] - optmodel.b[3] * bpChoices[optBP], optmodel.b[1], optmodel.b[2] + optmodel.b[3]]), 'exp_inoutflow_model2A.csv')
    return results, y_hat


    results = np.zeros((len(bpChoices)-1, 2))
    print bpChoices
    
    for i in range(len(bpChoices)-1):
        print i
        x2star = (predictor - bpChoices[i]) * np.greater(predictor, bpChoices[i])
        x1star = x1 * np.greater(predictor, bpChoices[i]) 
        tempPredictor = np.array(zip(x1, x1star, predictor, x2star))
        #fileLoc = filePath + 'temp.csv'
        #np.savetxt(fileLoc, tempPredictor, delimiter=',', fmt = '%s')
        #print tempPredictor
        tempmodel = ols.ols(response, tempPredictor,'y' ,predictorName)
        results[i,0] = i
        #results[i,1] = tempmodel.sse
        results[i,1] = tempmodel.R2

    optBP = int(results[np.argmax(results, axis = 0)[1],0])
    print 'Optimal Index:', optBP
    print 'Optimal changepoint: ', bpChoices[optBP], ' exp value: ', np.exp(bpChoices[optBP]), ' with R2 = ', results[optBP, 1]

    #x2star = (predictor - bpChoices[optBP]) * np.greater(predictor, bpChoices[optBP])
    #optPredictor = np.array(zip(predictor, x2star))
    #optmodel = ols.ols(response, optPredictor,'y',['x1', 'x2'])
    x1star = x1 * np.greater(predictor, bpChoices[optBP])
    x2star = (predictor - bpChoices[optBP]) * np.greater(predictor, bpChoices[optBP])
    optPredictor = np.array(zip(x1, x1star, predictor, x2star))
    optmodel = ols.ols(response, optPredictor,'y',predictorName)
    
    #return bpChoices[optBP], results, optmodel, optmodel.b[0]+optmodel.b[1]*predictor+optmodel.b[2]*x2star
    print results, optmodel.b
    print optmodel.summary()
    return results

def separateData(data1, data2, dataC, threshold):
    if (len(data1) <> len(data2)) or (len(data1) <> len(dataC)):
        print 'Error in input data length! --- seperateData.error.1'

    data1Left = []
    data1Right = []
    data2Left = []
    data2Right = []
    dataleftIdx = []
    datarightIdx = []
    
    for i in range(len(data1)):
        if dataC[i] < threshold:
            data1Left.append(data1[i])
            data2Left.append(data2[i])
            dataleftIdx.append(i)
        else:
            data1Right.append(data1[i])
            data2Right.append(data2[i])
            datarightIdx.append(i)

    return np.array(data1Left), np.array(data1Right), np.array(data2Left), np.array(data2Right), dataleftIdx, datarightIdx

def calculateR2(tleftmodel, trightmodel, totalmean):
    #self.sse = dot(self.e,self.e)/self.df_e
    #tsse = np.dot(tmodel.e, tmodel.e)/tmodel.df_e
    #tsse = np.dot(tmodel.e, tmodel.e)
    tsse = tleftmodel.sse * tleftmodel.df_e + trightmodel.sse * trightmodel.df_e
    
    #y_hat = np.dot(tmodel.x,tmodel.b)
    tsst = np.dot(tleftmodel.y - totalmean, tleftmodel.y - totalmean)
    tsst = tsst + np.dot(trightmodel.y - totalmean, trightmodel.y - totalmean)
    
    #print tleftmodel.sse * tleftmodel.df_e, trightmodel.sse * trightmodel.df_e, np.dot(tleftmodel.y - totalmean, tleftmodel.y - totalmean), np.dot(trightmodel.y - totalmean, trightmodel.y - totalmean)
    #print 1 - (tleftmodel.sse * tleftmodel.df_e)/np.dot(tleftmodel.y - totalmean, tleftmodel.y - totalmean)
    #print 1 - (trightmodel.sse * trightmodel.df_e)/np.dot(trightmodel.y - totalmean, trightmodel.y - totalmean)
    
    #return 1-tsse/tsst
    #print ((tleftmodel.nobs+trightmodel.nobs-1.0), (tleftmodel.nobs+trightmodel.nobs-tleftmodel.ncoef-trightmodel.ncoef)), ((tleftmodel.nobs+trightmodel.nobs-1.0)/(tleftmodel.nobs+trightmodel.nobs-tleftmodel.ncoef-trightmodel.ncoef))
    return 1 - (tsse/tsst)*((tleftmodel.nobs+trightmodel.nobs-1.0)/(tleftmodel.nobs+trightmodel.nobs-tleftmodel.ncoef-trightmodel.ncoef))

def pctmis(fobs, fpre):
    # Percentage Misallocated (Thomas, 1977)
    #print sum(abs(fobs-fpre)), sum(fobs)
    return 50*sum(abs(fobs-fpre))/sum(fobs)

def pctmis_whole(left_obs, left_pre, right_obs, right_pre):
    tsum = sum(abs(left_obs-left_pre)) + sum(abs(right_obs-right_pre)) + 0.0
    #print tsum, (sum(left_obs)+sum(right_obs))
    return 50*tsum/(sum(left_obs)+sum(right_obs))

def SRMSE(fobs, fpre):
    # Standardized Root Mean Square Error (Knudsen and Fotheringham, 1986)
    fe = fobs-fpre
    t = np.dot(fe, fe)/len(fobs)
    #print t
    t = t**0.5
    #print 't = ', t
    #print '(sum(fobs)/len(fobs)) = ', (sum(fobs)/len(fobs))
    return t/((sum(fobs)+0.0)/len(fobs))

def SRMSE_whole(left_obs, left_pre, right_obs, right_pre):
    lefte = left_obs - left_pre
    righte = right_obs - right_pre
    t = (np.dot(lefte, lefte) + np.dot(righte, righte))/(len(left_obs) + len(right_obs)) + 0.0
    t = t**0.5
    #print 't = ', t
    #print '(sum(fobs)/len(fobs)) = ', ((sum(left_obs) + sum(right_obs))/(len(left_obs)+ len(right_obs)))
    return t/((sum(left_obs) + sum(right_obs) + 0.0)/(len(left_obs)+ len(right_obs)))

def simpleRegression(response, predictor, predictorName):
    model = ols.ols(response, predictor,'y',predictorName)
    #model = ols.ols(response, predictor)
    print model.summary()
    print model.b
    print predictor
    #print 'b0 =',model.b[0], ', b1 =', model.b[1], ', bd =', model.b[2], ', MSE =', model.sse
    return np.dot(model.x, model.b), model

def separateRegression(response, predictor, sepData, bpChoices, predictorName, equaltestid, equaltestid2):
    results = np.zeros((len(bpChoices)-1, 4))
    print bpChoices
    chosencriterion = 'r2'  # max(r2) is the same criterion as min(sse)
    for bpid in range(len(bpChoices)-1):
        #print bpid
        responseLeft, responseRight, predictorLeft, predictorRight, dataleftIdx, datarightIdx = separateData(response, predictor, sepData, bpChoices[bpid])
        #print predictorLeft
        leftmodel = ols.ols(responseLeft, predictorLeft,'y',predictorName)
        rightmodel = ols.ols(responseRight, predictorRight,'y',predictorName)
        results[bpid,0] = bpid
        if chosencriterion == 'r2':           
            results[bpid,1] = leftmodel.R2adj
            results[bpid,2] = rightmodel.R2adj
            #results[bpid,3] = 1 - (leftmodel.e.var() + rightmodel.e.var())/(leftmodel.y.var() + rightmodel.y.var())
            results[bpid,3] = calculateR2(leftmodel, rightmodel, np.mean(response))
            #results[bpid,3] = (leftmodel.R2 + rightmodel.R2)/2
        elif chosencriterion == 'sse':
            results[bpid,1] = leftmodel.sse * leftmodel.df_e
            results[bpid,2] = rightmodel.sse * rightmodel.df_e
            results[bpid,3] = results[bpid,1] + results[bpid,2]
            #yhatL = np.dot(leftmodel.x, leftmodel.b)
            #print results[bpid,1], np.sum((responseLeft - yhatL)**2)

    print results
    if chosencriterion == 'r2':
        optBP = int(results[np.argmax(results, axis = 0)[-1],0])
    elif chosencriterion == 'sse':
        optBP = int(results[np.argmin(results, axis = 0)[-1],0])

    responseLeft, responseRight, predictorLeft, predictorRight, dataleftIdx, datarightIdx = separateData(response, predictor, sepData, bpChoices[optBP])
    leftmodel = ols.ols(responseLeft, predictorLeft,'y',predictorName)
    rightmodel = ols.ols(responseRight, predictorRight,'y',predictorName)

    #equaltestid = int(equaltestid)
    if equaltestid[0] > -1:
        temppredictorLeft = predictorLeft
        temppredictorLeft[:,equaltestid[1]] = temppredictorLeft[:,equaltestid[0]] + temppredictorLeft[:,equaltestid[1]]
        temppredictorLeft = np.delete(temppredictorLeft, equaltestid[0], 1)
        temppredictorRight = predictorRight
        temppredictorRight[:,equaltestid[1]] = temppredictorRight[:,equaltestid[0]] + temppredictorRight[:,equaltestid[1]]
        temppredictorRight = np.delete(temppredictorRight, equaltestid[0], 1)
        temppredictorName = np.delete(predictorName, equaltestid[0], None)
        templeftmodel = ols.ols(responseLeft, temppredictorLeft,'y',temppredictorName)
        temprightmodel = ols.ols(responseRight, temppredictorRight,'y',temppredictorName)
        fleft = (leftmodel.R2 - templeftmodel.R2) * (leftmodel.nobs - len(predictorName) - 2) / (1 - leftmodel.R2)
        fright = (rightmodel.R2 - temprightmodel.R2) * (rightmodel.nobs - len(predictorName) - 2) / (1 - rightmodel.R2)
        pleft = 1 - f.cdf(fleft, 1, leftmodel.nobs - len(predictorName) - 2)
        pright = 1 - f.cdf(fright, 1, rightmodel.nobs - len(predictorName) - 2)

        if equaltestid2[0] > -1:
            temppredictorLeft = predictorLeft
            temppredictorLeft[:,equaltestid2[1]] = temppredictorLeft[:,equaltestid2[0]] + temppredictorLeft[:,equaltestid2[1]]
            temppredictorLeft = np.delete(temppredictorLeft, equaltestid2[0], 1)
            temppredictorRight = predictorRight
            temppredictorRight[:,equaltestid2[1]] = temppredictorRight[:,equaltestid2[0]] + temppredictorRight[:,equaltestid2[1]]
            temppredictorRight = np.delete(temppredictorRight, equaltestid2[0], 1)
            temppredictorName = np.delete(predictorName, equaltestid2[0], None)
            templeftmodel = ols.ols(responseLeft, temppredictorLeft,'y',temppredictorName)
            temprightmodel = ols.ols(responseRight, temppredictorRight,'y',temppredictorName)
            fleft = (leftmodel.R2 - templeftmodel.R2) * (leftmodel.nobs - len(predictorName) - 2) / (1 - leftmodel.R2)
            fright = (rightmodel.R2 - temprightmodel.R2) * (rightmodel.nobs - len(predictorName) - 2) / (1 - rightmodel.R2)
            pleft2 = 1 - f.cdf(fleft, 1, leftmodel.nobs - len(predictorName) - 2)
            pright2 = 1 - f.cdf(fright, 1, rightmodel.nobs - len(predictorName) - 2)

    yhatL = np.dot(leftmodel.x, leftmodel.b)
    yhatR = np.dot(rightmodel.x, rightmodel.b)
    yhat = np.zeros(len(response))

    for i in range(len(yhatL)):
        yhat[dataleftIdx[i]] = yhatL[i]

    for i in range(len(yhatR)):
        yhat[datarightIdx[i]] = yhatR[i]

    yhat = np.exp(yhat)
    fileLoc = filepath + 'separateR_model2_y_hat.csv'
    #np.savetxt(fileLoc, yhat, delimiter=',', fmt = '%s')
    #print 'Optimal Index:', optBP
    #print 'Optimal changepoint: ', bpChoices[optBP], ' exp value: ', np.exp(bpChoices[optBP]), ' with R2 = ', calculateR2(leftmodel, rightmodel, np.mean(response))
   
    print '----------------------------- left model -----------------------------'
    print leftmodel.summary()
    print '----------------------------- right model -----------------------------'
    print rightmodel.summary()

    print 'Optimal Index:', optBP
    print 'Optimal changepoint: ', bpChoices[optBP], ' exp value: ', np.exp(bpChoices[optBP]), ' with adjR2 = ', results[optBP, -1]
   
    outputstring = 'before bp'
    for i in range(len(predictorName)+1):
        outputstring += ', b' + str(i) + ' = ' + "%.2f" %(leftmodel.b[i]) + '(' + "%.3f" %(leftmodel.se[i])+ ')'
    outputstring += ', with adjR2 = ' + "%.4f" %(leftmodel.R2adj)
    if equaltestid[0] > -1:
        outputstring += ', b1 <> b2 with pvalue = ' + "%.3f" %(pleft)
        if equaltestid2[0] > -1:
            outputstring += ', b3 <> b4 with pvalue = ' + "%.3f" %(pleft2)
    print outputstring

    outputstring = 'after bp'
    for i in range(len(predictorName)+1):
        outputstring += ', b' + str(i) + ' = ' + "%.2f" %(rightmodel.b[i]) + '(' + "%.3f" %(rightmodel.se[i])+ ')'
    outputstring += ', with adjR2 = ' + "%.4f" %(rightmodel.R2adj)
    if equaltestid[0] > -1:
        outputstring += ', b1 <> b2 with pvalue = ' + "%.3f" %(pright)
        if equaltestid2[0] > -1:
            outputstring += ', b3 <> b4 with pvalue = ' + "%.3f" %(pright2)
    print outputstring

    syhat, smodel = simpleRegression(response, predictor, predictorName)

    fsmr = (results[optBP, -1] - smodel.R2) * (smodel.nobs - len(predictorName) - 2) / (1 - results[optBP, -1])
    psmr = 1 - f.cdf(fsmr, 1, smodel.nobs - len(predictorName) - 2)

    print '--------------------Simple Regression Result-----------------------'
    outputstring = ''
    for i in range(len(predictorName)+1):
        outputstring += 'b' + str(i) + ' = ' + "%.2f" %(smodel.b[i]) + '(' + "%.3f" %(smodel.se[i])+ '), '
    outputstring += 'with adjR2 :' + '%.4f' %(smodel.R2adj) + ', separate regression with pvalue = ' + "%.3f" %(psmr)
    print outputstring

    #calpredictedvalue(predictor, bpChoices[optBP], zip(leftmodel.b, rightmodel.b), 'exp_inoutflow_model2B.csv')
    #calconfidenceinterval(predictorLeft, predictorRight, [leftmodel.sse, rightmodel.sse], response, predictor, bpChoices[optBP], zip(leftmodel.b, rightmodel.b), 'ci_model2B.csv')
    return results, yhat

def separatePoisson(response, predictor, sepData, bpChoices):
    results = np.zeros((len(bpChoices)-1, 3))
    print bpChoices
    chosencriterion = 'PM'  # PM ; SRMSE
    for bpid in range(len(bpChoices)-1):
        #print bpid
        responseLeft, responseRight, predictorLeft, predictorRight, dataleftIdx, datarightIdx = separateData(response, predictor, sepData, bpChoices[bpid])
        #print predictorLeft
        #inputPredictor = sm.add_constant(inputPredictor, prepend=False)
        leftmodel  = sm.Poisson(responseLeft, predictorLeft)        
        leftmodel.solution = leftmodel.fit(method="bfgs")
        
        rightmodel = sm.Poisson(responseRight, predictorRight)
        rightmodel.solution = rightmodel.fit(method="bfgs")
        
        #leftmodel = ols.ols(responseLeft, predictorLeft,'y',predictorName)
        #rightmodel = ols.ols(responseRight, predictorRight,'y',predictorName)
        leftpre = leftmodel.solution.predict()
        rightpre = rightmodel.solution.predict()
        
        results[bpid,0] = bpid          
        results[bpid,1] = pctmis_whole(responseLeft, leftpre, responseRight, rightpre) # PM
        results[bpid,2] = SRMSE_whole(responseLeft, leftpre, responseRight, rightpre) # SRMSE

    print results
    if chosencriterion == 'PM': #lower better
        optBP = int(results[np.argmin(results, axis = 0)[1],0])
    elif chosencriterion == 'SRMSE': # lower better
        optBP = int(results[np.argmin(results, axis = 0)[2],0])

    responseLeft, responseRight, predictorLeft, predictorRight, dataleftIdx, datarightIdx = separateData(response, predictor, sepData, bpChoices[optBP])
    leftmodel  = sm.Poisson(responseLeft, predictorLeft)        
    leftmodel.solution = leftmodel.fit(method=poissonmethod)
        
    rightmodel = sm.Poisson(responseRight, predictorRight)
    rightmodel.solution = rightmodel.fit(method=poissonmethod)

    yhatL = leftmodel.solution.predict()
    yhatR = rightmodel.solution.predict()
    yhat = np.zeros(len(response))

    for i in range(len(yhatL)):
        yhat[dataleftIdx[i]] = yhatL[i]

    for i in range(len(yhatR)):
        yhat[datarightIdx[i]] = yhatR[i]

    yhat = np.exp(yhat)
    fileLoc = filepath + 'separateR_model2_y_hat.csv'
    #np.savetxt(fileLoc, yhat, delimiter=',', fmt = '%s')
    #print 'Optimal Index:', optBP
    #print 'Optimal changepoint: ', bpChoices[optBP], ' exp value: ', np.exp(bpChoices[optBP]), ' with R2 = ', calculateR2(leftmodel, rightmodel, np.mean(response))
   
    print '----------------------------- left model -----------------------------'
    print leftmodel.solution.summary()
    print '----------------------------- right model -----------------------------'
    print rightmodel.solution.summary()

    print 'Optimal Index:', optBP
    print 'Optimal changepoint: ', bpChoices[optBP], ' with PM = ', results[optBP,1], ' and SRMSE = ', results[optBP,2]

    '''
    outputstring = 'before bp'
    for i in range(len(predictorName)+1):
        outputstring += ', b' + str(i) + ' = ' + "%.2f" %(leftmodel.b[i]) + '(' + "%.3f" %(leftmodel.se[i])+ ')'
    outputstring += ', with adjR2 = ' + "%.4f" %(leftmodel.R2adj)
    if equaltestid[0] > -1:
        outputstring += ', b1 <> b2 with pvalue = ' + "%.3f" %(pleft)
        if equaltestid2[0] > -1:
            outputstring += ', b3 <> b4 with pvalue = ' + "%.3f" %(pleft2)
    print outputstring

    outputstring = 'after bp'
    for i in range(len(predictorName)+1):
        outputstring += ', b' + str(i) + ' = ' + "%.2f" %(rightmodel.b[i]) + '(' + "%.3f" %(rightmodel.se[i])+ ')'
    outputstring += ', with adjR2 = ' + "%.4f" %(rightmodel.R2adj)
    if equaltestid[0] > -1:
        outputstring += ', b1 <> b2 with pvalue = ' + "%.3f" %(pright)
        if equaltestid2[0] > -1:
            outputstring += ', b3 <> b4 with pvalue = ' + "%.3f" %(pright2)
    print outputstring

    syhat, smodel = simpleRegression(response, predictor, predictorName)

    fsmr = (results[optBP, -1] - smodel.R2) * (smodel.nobs - len(predictorName) - 2) / (1 - results[optBP, -1])
    psmr = 1 - f.cdf(fsmr, 1, smodel.nobs - len(predictorName) - 2)

    print '--------------------Simple Regression Result-----------------------'
    outputstring = ''
    for i in range(len(predictorName)+1):
        outputstring += 'b' + str(i) + ' = ' + "%.2f" %(smodel.b[i]) + '(' + "%.3f" %(smodel.se[i])+ '), '
    outputstring += 'with adjR2 :' + '%.4f' %(smodel.R2adj) + ', separate regression with pvalue = ' + "%.3f" %(psmr)
    print outputstring
    '''
    #calpredictedvalue(predictor, bpChoices[optBP], zip(leftmodel.b, rightmodel.b), 'exp_inoutflow_model2B.csv')
    #calconfidenceinterval(predictorLeft, predictorRight, [leftmodel.sse, rightmodel.sse], response, predictor, bpChoices[optBP], zip(leftmodel.b, rightmodel.b), 'ci_model2B.csv')
    return results, yhat

def calconfidenceinterval(leftx, rightx, mse, y, x, bp, b,filename):
    leftx = np.array(zip(np.ones(len(leftx)), leftx[:,0], leftx[:,1]))
    rightx = np.array(zip(np.ones(len(rightx)), rightx[:,0], rightx[:,1]))

    templeftx = np.linalg.inv(np.dot(leftx.T, leftx))
    temprightx = np.linalg.inv(np.dot(rightx.T, rightx))

    b = np.array(b)
    t = [stats.t.ppf(0.975, len(leftx) - len(leftx[0,:])), stats.t.ppf(0.975, len(rightx) - len(rightx[0,:]))]

    yhat = []
    i = 0
    #print len(x)
    for item in x:
        if item[1] < bp:
            bpflag = 0
            tempx = templeftx
        else:
            bpflag = 1
            tempx = temprightx
        temp = b[0, bpflag] + b[1, bpflag] * item[0] + b[2, bpflag] * item[1]
        xh = np.array([1, item[0], item[1]])
        #print np.dot(np.dot(xh, tempx), xh.T), mse[bpflag]*(np.dot(np.dot(xh, tempx), xh.T)
        se = math.pow(mse[bpflag]*(np.dot(np.dot(xh, tempx), xh.T)), 0.5) * t[bpflag]
        pise = math.pow(mse[bpflag]*(1 + np.dot(np.dot(xh, tempx), xh.T)), 0.5) * t[bpflag]
        cil = 0
        cih = 0
        pil = 0
        pih = 0
        if y[i] < (temp-se):
            cil = 1
        if y[i] > (temp+se):
            cih = 1
        if y[i] < (temp-pise):
            pil = 1
        if y[i] > (temp+pise):
            pih = 1
        yhat.append([temp, temp-se, temp+se, temp-pise, temp+pise,cil,cih,pil,pih])
        i += 1

    yhat = np.exp(yhat)
    fileLoc = filepath + filename   #'exp_inoutflow_model2A.csv'
    headerstr = 'yhat,cileft,ciright,pileft,piright,cil,cih,pil,pih'
    #print fileLoc
    np.savetxt(fileLoc, yhat, delimiter=',', header = headerstr, fmt = '%f')
    
    #plt.plot(range(len(y)), y, 'g', label='Observed')
    #plt.plot(range(len(y)), yhat[:,0], 'b', label='Predicted')
    #plt.plot(range(len(y)), yhat[:,1], 'b', label='CI Left')
    #plt.plot(range(len(y)), yhat[:,2], 'b', label='CI Right')
    #plt.plot(range(len(y)), yhat[:,3], 'r', label='PI Left')
    #plt.plot(range(len(y)), yhat[:,4], 'r', label='PI Right')
    #plt.show()
    

def calpredictedvalue(x, bp, b, filename):
    # x: [f1f2, dist]
    # b = [[b00, b10], [b01, b11], [b02, b12]]
    print b
    yhat = []
    b = np.array(b)
    bpflag = 1  # 0: before, 1: after
    for item in x:
        if item[1] < bp:
            bpflag = 0
        else:
            bpflag = 1
        temp = b[0, bpflag] + b[1, bpflag] * item[0] + b[2, bpflag] * item[1]
        #temp = np.dot(item, b[:,bpflag])
        yhat.append(temp)

    yhat = np.array(yhat)
    yhat = np.exp(yhat)
    #print yhat
    #fileLoc = filepath + 'model2b_yhat.csv'
    #np.savetxt(fileLoc, yhat, delimiter=',', fmt = '%f')
    expinoutflow = np.zeros((len(inoutflow), 2)) # [inflow, outflow, obs-exp inflow, obs-exp outflow]
    #print len(yhat), len(data)
    #print zip(data, yhat)
    i = 0
    for item in data:
        expinoutflow[int(item[1]), 1] += yhat[i]
        expinoutflow[int(item[2]), 0] += yhat[i]
        i += 1

    fileLoc = filepath + filename   #'exp_inoutflow_model2A.csv'
    headerstr = 'inflow,outflow'
    #np.savetxt(fileLoc, expinoutflow, delimiter=',', header = headerstr, fmt = '%f')
    
def dataprepare():
    filepath = 'C:/_DATA/migration_census_2000/REDCAPLUS_100Region/'

    #dataCSV = filepath + 'migration_100regs_dist.csv'
    dataCSV = filepath + 'census_migration_100regs_dis.csv'
    data = build_data_list(dataCSV)

    print data

    file = filepath + 'census_migration_100regs_pop.csv'
    #file = filepath + '100reg_pop.csv'
    pop = build_data_list(file)

    file = filepath + '100reg_popdensity.csv'
    tdata = build_data_list(file, "float")
    popden = tdata[:,-1]

    '''
    grossflowCSV = filepath + 'migration_100regs_inoutflow.csv'
    #grossflowCSV = filepath + 'census_migration_100regs_inoutflow.csv'
    grossflow = build_data_list(grossflowCSV)

    inflow = np.zeros(len(grossflow))
    outflow = np.zeros(len(grossflow))

    for item in data:
        outflow[item[1]] += item[-1]
        inflow[item[2]] += item[-1]

    print zip(inflow, outflow)
    '''
        
    output = []
    
    for item in data:
        #dist = int(cal_dist(centroid[item[0], 0], centroid[item[0], 1], centroid[item[1], 0], centroid[item[1], 1]))
        #output.append([item[0], grossflow[item[1], 1], grossflow[item[2], 0], pop[item[1], 1], pop[item[2], 0], item[3]])
        output.append([item[0]/1000, pop[item[1], 1], pop[item[2], 0], popden[item[1]], popden[item[2]], item[3]])

    #print output
    headerstr = 'dist,p1,p2,pd1,pd2,flowVol'
    print headerstr
    np.savetxt(filepath + '100reg_regdata_popden.csv', output, delimiter=',', header = headerstr, fmt = '%s')
    return output
    
#--------------------------------------------------------------------------
#MAIN

if __name__ == "__main__":
    print '===================================================='
    print "begin at " + getCurTime()

    filepath = 'C:/_DATA/migration_census_2000/REDCAPLUS_100Region/'
    #filepath = 'C:/_DATA/migration_census_2000/state/'
    #file = filepath + 'state_regressiondata.csv'
    file = filepath + '100reg_regressiondata.csv' #[dist,f1,f2,p1,p2,vol]: 1:out, 2:in
    #file = filepath + 'regressiondataeast.csv'
    #file = filepath + 'regressiondatawest.csv'    
    data = build_data_list(file)   

    #data = dataprepare()    #[dist,p1,p2,pd1,pd2,flowVol]: 1:out, 2:in
    data = np.array(data)

    lndata = np.log(data)

    dist = data[:,0]
    lndist = lndata[:,0]

    inputPredictor = []
    #i = 0
    for i in range(len(data)):
        inputPredictor.append([lndata[i,3], lndata[i,4], lndata[i,0]])
    inputPredictor = np.array(inputPredictor)

    #rand_data = sm.datasets.randhie.load()
    print 'Poisson Regression Started...'
    inputPredictor = sm.add_constant(inputPredictor, prepend=False)
    poisson_mod = sm.Poisson(data[:,-1], inputPredictor)
    poisson_res = poisson_mod.fit(method="powell")
    print poisson_res.summary()
    print 'PM = ', pctmis(data[:,-1], poisson_res.predict())
    print 'SRMSE = ', SRMSE(data[:,-1], poisson_res.predict())
    
    #inputPredictor = lndata[:,[3,4,0]]

    bpChoices = range(int(min(dist)/10+20)*10,int(max(dist)/10-100)*10,10)
    #bpChoices = np.log(range(90,4200,10))
    #bpChoices = np.log(range(90,2100,10))  #east
    #bpChoices = range(90,2100,10)  #east
    #bpChoices = np.log(range(140,3100,10))  #west
    #bpChoices = range(140,3100,10)  #west

    inputPredictorName = ['lnp1', 'lnp2', 'dist']
    poissonmethod = 'bfgs' # 'newton' 'nm' ' bfgs' 'powell' 'cg' 'ncg'
    
    '''
    responseLeft, responseRight, predictorLeft, predictorRight, dataleftIdx, datarightIdx = separateData(data[:,-1], inputPredictor, dist, 890)
    leftmodel  = sm.Poisson(responseLeft, predictorLeft)        
    leftmodel.solution = leftmodel.fit(method=poissonmethod)
        
    rightmodel = sm.Poisson(responseRight, predictorRight)
    rightmodel.solution = rightmodel.fit(method=poissonmethod)
    
    leftpre = leftmodel.solution.predict()
    rightpre = rightmodel.solution.predict()

    ttt1 = pctmis_whole(responseLeft, leftpre, responseRight, rightpre) # PM
    ttt2 = SRMSE_whole(responseLeft, leftpre, responseRight, rightpre) # SRMSE

    print '----------------------------- left model -----------------------------'
    print leftmodel.solution.summary()
    print '----------------------------- right model -----------------------------'
    print rightmodel.solution.summary()

    print 'PM = ', ttt1
    print 'SRMSE = ', ttt2
    '''
    

    #inputPredictorName = ['lnp1p2', 'lndist']
    #yhat, smodel = simpleRegression(lndata[:,-1], inputPredictor, inputPredictorName)
    #inputPredictorName = ['lnf1', 'lnf2', 'lnf12', 'lnf22', 'lndist']
    #result, yhat = separateRegression(lndata[:,-1], inputPredictor, dist, bpChoices, inputPredictorName, [0,1], [-1,3])
    #result, yhat = separatePoisson(data[:,-1], inputPredictor, dist, bpChoices)
    #inputPredictorName = ['lnp1', 'lnp2', 'lndist']
    #results, y_hat = pickBreakpoint(lndata[:,-1], inputPredictor, lndist, bpChoices, inputPredictorName)

    #fileLoc = filepath + 'piecewiseR_model2_results.csv'
    #fileLoc = filepath + 'separateR_model2_results.csv'
    headerstr = 'idx,PM,RSMSE'
    fileLoc = filepath + 'temp.csv'
    #np.savetxt(fileLoc, result, delimiter=',', header = headerstr, fmt = '%s')
    #np.savetxt(fileLoc, poisson_res.predict(), delimiter=',', header = headerstr, fmt = '%s')
    